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Abstract 

We present a set of molecular dynamics (MD) simulations of strongly 
charged, flexible polyelectrolyte chains under poor solvent conditions 
in a salt free solution. Structural properties of the chains and of the 
solutions are reported. By varying the polymer density and the electro- 
static interaction strength we study the crossover from a dominating 
electrostatic interaction to the regime of strong screening, where the 
hydrophobic interactions dominate. During the crossover a multitude 
of structures is observed. In the limit of low polymer density strongly 
stretched, necklace like conformations are found. In the opposite limit 
of high polymer density which is equivalent to strongly screened elec- 
trostatic interactions, we find that the chains are extremely collapsed, 
however we observe no agglomeration or phase separation. The inves- 
tigations show that the density of free charges is one of the relevant 
parameters which rules the behavior of the system and hence should 
be used as a parameter to explain experimental results. 
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1 Introduction 



Unlike neutral polymers ^, [| , the understanding of the behavior of elec- 
trically charged macromolecules, short polyelectrolytes, is still rather poor. 
The long range nature of the electrostatic interactions introduces new length 
and time scales that render the analytical description very complicated [Q] 
and prevent simple scaling arguments to work equally successfully as in the 
neutral case. Usually, analytic theory || assumes ^-properties for the un- 
charged monomers, while simulations use either 9- or good solvent conditions 
|], [/], ||, ^, 10 1 . Nevertheless, most experimentally relevant systems have or- 



ganic, non-polar backbones in polar solvents (water) so that poor solvent 
conditions are more appropriate [12]. The analytic description of these sys- 
tems is not far developed. In a first attempt, Dobrynin et al. [ 13 1 studied 
a single weakly charged chain under poor solvent conditions. Starting from 
the description of a charged drop, originally given by Lord Rayleigh in 1882 
[14 1, their theory predicts that when the charge density or the strength of 
the electrostatic repulsion increases, the globular state becomes unstable 
and starts to split into smaller globules, which try to separate. Due to the 
covalent bonds the parts cannot dissociate completely. By maximizing their 
distance, they form structures which are called pearl-necklace chains. 

The approach of Dobrynin et al. is supported by a Monte Carlo simula- 
tion of a single chain that shows a cascade from one to two and three globules 
with increasing strength of the electrostatic repulsion In their study 

every monomer carries a fractional charge, the polyelectrolyte is weakly 
charged and no explicit counterions are taken into account. Additionally 
to this previous Monte Carlo study, there exist some experimental evidence 
[16 1 supporting the picture of Dobrynin at al. A recent Monte-Carlo simula- 
tion of single Debye-Hiickel chains near the 6 transition also supports these 
findings [fT^| . However, there has been also put some doubt on the validity 
of this scheme in the case of discrete ions p8| . Our study will use many 
chains with discrete counterions, where the charges all interact via the full 
Coulomb potential, simulated wit molecular dynamics (MD). 

Alternatively, a cylindrical instability was proposed by Khokhlov and his 
collaborators (i9[| . They argue that there will be a micro-phase separation 
due to a kind of agglomeration of the collapsed chains. This state would 
increase the entropy of the condensed counterions which could explore the 
whole volume of the phase separated chains. In an earlier short communi- 
cation |^0| we showed that based on simulational data no phase separation 
was observed and gave a simple argument for a stable colloid like state of 
single collapsed chains. In this article we will present in detail our simula- 
tions of polyelectrolytes in poor solvents. The rest of the paper is organized 
as follows: In section ||| we will describe our model, and show in section ^ 
the strong influence of the hydrophobic interaction, when it is put into a 
right proportion to the electrostatic interaction. In section || we will present 
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the single chain properties of a polyelectrolyte under poor solvent condition, 
which shows indeed a cascade of pearl-necklace structures, if one varies the 
polymer density. We also explain why the colloid like state of the collapsed 
chain seems to be the stable one, and why no phase separation is observed. 
As last point in section [5] we compare our results to experimental data and 
suggest for future experiments those parameters which we consider to be the 
important quantities to vary. 



2 Model 

The polyelectrolyte chains are modeled as bead spring chains of Lennard- 
Jones (LJ) particles. 

For good solvent conditions, a shifted LJ potential is used to describe 
the purely repulsive excluded volume interaction between all monomers: 



ULJ,hydrophil(r) 



eLj[(f) 12 - (f) 6 ] + e LJ for r < r min = 2§c7 
for r > r min = 2$o 



The parameter a is used to set the length scale. The energy parameter 
€lj controls the strength of the interaction, and its value for the good sol- 
vent case is fixed to €lj = 1.0/^bT, where ks is the Boltzmann constant 
and T is the temperature. Polymers that are in good solvent and whose 
monomers interact via this potential are called "hydrophilic" for the rest of 
this paper. Although we do not have any explicit solvent and no specific 
potentials for the solvent water, we use this abbreviation because water is in 
experiments, technical applications, and nature the most common solvent 
for polyelectrolytes. As a consequence, particles under poor solvent con- 
ditions are called "hydrophobic". Because we do not consider any explicit 
interactions between solutes and solvent molecules one way to model poor 
solvent conditions is through an additional attractive interaction between 
the chain monomers. 

We represent this potential by a standard LJ potential with a larger 
cut-off (R c = 2.5a): 

J* ■ e L j[(f) 12 - (f ) 6 - c(R c )] for r < R c 

U LJ,hydrophob\ r ) — \ (2) 
[0 for r > R c 

The function c(R c ) is chosen as c(R c ) = {-§£} 12 — (j|:) 6 to give a potential 
value of zero at the cutoff. This potential is used only for chain beads, 
because we assume that the counterions do not have any hydrophobic parts 
in their interaction, as is reasonable for alkali metals, and describe their 
interaction by equation |l|. 
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The connectivity of the bonded monomers is assured by the FENE (finite 
extension nonlinear elastic) potential: 

1 r 2 
C/feneO) = --fcFENE^oM 1 ~~ 7^) ( 3 ) 

where r is the distance of the interacting particles. The parameters are set 
the following way: the spring constant /cfene = 7.0^- and i?o = 2<r. Rq 
gives the maximum extension of the bond, at which the interaction energy 
becomes infinite. This potential gives rather stiff bond lengths with fluc- 
tuations being smaller than 10 %. The nonlinearity allows a very efficient 
mixing of the different vibrational modes. 

The Bjerrum length is a measure of the strength of the electrostatic 
interaction, which is defined as the length at which the electrostatic energy 
equals the thermal energy: 

^ B 47re 5 e A:BT' ^ 

where e is the unit charge of the interacting particles, and cq and 65 are the 
permittivity of the vacuum and of the solvent, respectively (As = 7.1 A in 
water). Here we restrict ourselves to monovalent species in salt free solu- 
tion. The solvent is taken into account via the dielectric background, whose 
properties influence the form (hydrophil or hydrophob) of the interaction 
potential of the monomers. Unless mentioned the Bjerrum length is always 
fixed to Xb = 3.0a. 

No assumptions concerning screening effects are made, instead we are 
using the pure Coulomb energy to calculate the interactions between all 
charged particles: 

E^{r ij ) = k B T ) ^i, (5) 

where qi = +1 for the charged chain monomers, and qi = — 1 for the coun- 
terions, which are treated explicitly as charged LJ particles. 

The central simulation box contained 16 or 32 chains subject to periodic 
boundary conditions. Each chain consisted of 94 monomers. The fraction 
/ of charged monomers on the chain is / = |, meaning that every third 
monomer carries a unit charge. We also put a charge on the chain ends, 
hence the number of charges per chain is always 32. A velocity Verlet algo- 
rithm is used to integrate the equation of motion. The system is coupled to 
a thermostat via a standard Langevin equation with a friction term and a 
stochastic force. The friction coefficient V is set to Y = t~ 1 and the time step 
At varies between At = 0.003r and 0.0125r depending on e^j, where r is 
the standard L J time unit |7| . Because of the periodic boundary conditions 
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the Coulomb interaction is calculated with a fast variant of the Ewald sum- 



mation [21], the Particle-Mesh-Ewald (PME) algorithm i n the version 
presented by Petersen 23]. The basic idea is to calculate the Fourier part 
of an Ewald summation with a Fast Fourier Transform on a mesh to speed 
up the calculation. Our version is more efficient than the highly optimized 
version of the Ewald sum of Ref.Q if more than 800 charges are simulated. 
The parameters of the PME method were always chosen such that the re- 
sulting errors in the electrostatic forces were well below the random forces 
produced by the Langevin heat bath. At the time of the production of our 
data this method was the only mesh based Ewald sum, for which analytic 
error estimates |2jJ were available. Note, however, that meanwhile a better 
method, also with analytic error estimates, is available 24]. 

Our polyelectrolyte systems are strongly charged. By this we mean that 
the Manning parameter which is defined as the Bjerrum length divided by 
the length per unit charge, £m = -gf, is close to the critical value. The 
critical Manning parameter, at which counterion condensation is commonly 
expected to set in, is equal to £^ = 1 for monovalent counterions. In 
our simulation, the average bond length is given by b ~ 1.08<t, leading 
to £jyf sa 0.927, so that we are just below ^* M , therefore the expression 
"strongly" charged. 

Our model polyelectrolyte system can be mapped to the "standard poly- 
electrolyte" sulfonated polystyrene (PSS), as can be seen for instance from 
the Manning parameter. A PSS-monomer extends over roughly 2.5A and 
the Bjerrum length in water is about 7.1 A leading to about one charge per 
three monomers. Each chain contains 94 monomers meaning 32 charges and 
models a molecular weight of roughly 20.000 -^-j, if a sodium counterion is 
assumed. This molecular weight, though not very large, is well within the 
experimentally analyzed regime |12| |. 

A crucial parameter in the simulations will be the concentration of par- 
ticles in the simulation box of length L. No difference will be made be- 
tween counterions, charged and uncharged monomers because they are all 
Lennard-Jones-particles. This concentration will be called density for the 
rest of this paper and therefore be denoted by p. By varying the box length 
L we can vary the density from 2 • lO^V" 3 to 2 • lO" 6 ^ 3 , which is cover- 
ing more than the whole experimentally accessible range from very dense to 
extremely dilute solutions. 



3 Influence of the hydrophobic interaction strength 



To characterize the hydrophobicity of a polyelectrolyte chain one has to 
develop a measure for the relative strengths of the hydrophobic and the 
electrostatic interaction. The strength of the electrostatic interaction is de- 
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termined by Xb- The attractive interaction in Lennard- Jones systems (as 
for instance neutral polymers) is usually tuned by variation of the temper- 
ature T. The aim of the present work is to study the influence of the poor 
solvent with fixed boundary conditions (constant Xb and constant T). For 
the case of polyelectrolytes a variation of the temperature in the simulations 
changes simultaneously the strength of the electrostatics for fixed Xb- For 
this reason, the hydrophobicity is varied by changing the parameter clj in 
the LJ potential (||) as an independent variable. We fix the energy scale by 
setting the thermal energy equal to ksT = 1.0. In a simulation this is a 
common way to separate and isolate effects. We are well aware of the fact 
that in experiments a change of temperature influences in a rather complex 
manner all relevant quantities. 

In order to make contact with neutral polymer theory the 9 point for the 
neutral polymer system has to be determined first. The 9 point is usually 
defined by the temperature at which the second virial coefficient between 
the chains vanishes. There the attractive and repulsive interactions cancel 
up to second order, so that the polymers show (almost) ideal random walk 
behavior. For the underlying neutral LJ system only a crude estimate for 
the 9 point has been found so far in Ref. p5[| . The temperature and the 
LJ energy parameter e£j of this point were determined as li > 3.0 ± 0.2 

(n) 

for € L j = 1.0. Converted to our notation this is equivalent to a value of 

t£je — 3 a ^ a thermal energy ksT = 1. To get a more precise estimate on 
the location of the 9 point we therefore employed a scaling approach that 
was shown to be very successful for polymers and star-polymers on a face- 
centered cubic lattice (fee) |^[ . This approach utilizes that the end-to-end 
distance can be expressed as: 

R 2 end = N-g(f -m), (6) 

where g is an unknown scaling function, and iV is the number of monomers 
of a chain. The argument T = T ^ Tg describes the reduced distance to the 

critical point. This relation is easily translated to an e^j-dependent formula. 

R 2 

The main point is that for fixed values of the ratio the argument of the 
scaling function must be also fixed. It follows that 

f ~ NT* (7) 

_1 B? 

Plotting T versus N 2 for fixed ratio of — sp* results in straight lines that 

should converge to the 9 point for ^ 0, neglecting some deviations due 
to statistical errors and higher order corrections to scaling. Neutral chains 
from N = 32 up to N = 2048 monomers were studied. The scaling plot can 
be viewed in Figure [l], showing that the effective interaction parameter for 
9 conditions is given by e^J e = [2.96 ± 0.06]" 1 = 0.34 ± 0.02. 
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Based on this information for the neutral system three simulations of the 
charged system were performed at a fixed density of p = 10 _4 cj~ 3 , but for 
the three different values e^j = 0.25, 0.5, and 1.0. While e^j = 0.25 is well 
above the neutral 6 point, where a behavior similar to the hydrophilic system 
can be expected, e^j = 0.5 should show some signs of the hydrophobic 
influence. For that reason, the system with e^j = 0.5 will be called "weakly 
hydrophobic" for the rest of the paper, while the system with eu = 1-0 is 
deep in the hydrophobic regime. 

We first looked at the spherically averaged structure factor of a single 
chain 



S sp (q) = ( ^7 



^exp(i<f • (fi -fj)) 



(8) 



and its component along the first principal axis of inertia, Sn, which is 
obtained by restricting the q - vectors to the direction of the calculated 
principal axis. The values for S sp (q) and Su(q) are displayed in the top and 
bottom part of Figure |2|, respectively, where it is found that all three of these 
systems are very similar and do not differ significantly from the known be- 
havior of hydrophilic polyelectrolytes. Because the charge density strongly 
influences the screening of the electrostatic interactions, the comparison of 
charged systems for different values of €lj is extended over a whole range of 
densities. In Figure || the radius of gyration Rq is displayed versus various 
densities. The attractive potentials cause locally more compact conforma- 
tions, leading to smaller Rq, which can also be seen in Figure |], where the 
spherically averaged structure factor is displayed for p = 10 _1 <7 -3 . 

The main results of these preliminary investigations show, that the elec- 
trostatic interaction is still too strong compared to the influence of the 
hydrophobicity. However, under certain conditions, even the "weak" hy- 
drophobicity has pronounced effects. 

One of theses conditions is the relative strength of the electrostatic inter- 
action. The strength of the Coulomb interaction can be conveniently altered 
by varying the Bjerrum length A#. It is well known that if one increases Xb 
by starting from small values, this will lead to a non-monotonic behavior of 



the chain dimensions [Q, which is due to counterion condensation [27, 28, 2£]. 



The basic idea of counterion condensation is a generalization of the well 



known solution of Onsager [30] for the infinitely long and thin line charge. 
If the energy density on the rod exceeds a certain value it becomes more 
favorable to catch one of the free ions to neutralize one of the charges on the 
chain. The calculation results in the condition that the distance between two 
adjacent unit charges on the chain has to be bigger than A#, which means 
that their mutual pure Coulomb repulsion energy has to be smaller than 
&bT. For that reason, the entropy of the counterion, that is lost in the con- 
densation process, is assumed to equal ksT. By analogy, this argument is 
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applied to any (flexible) polyelectrolyte, which is, however, questionable. We 
define a counterion as condensed, if its total energy is smaller than —lksT. 
This way is easy to handle and strongly supported by a visual inspection of 
the set of condensed counterions of our conformations. It has some advan- 
tages over the geometrical definition, where one looks at the distance of a 
counterion to its nearest polymer, and sets some arbitrary bound, because 
an ion can feel a strong attraction to several charged monomers even at a 
relatively large distance. 

The results of Ref. |?J and our new data show that first, increasing Xb 
stretches the chains through the increased repulsion of the equally charged 
chain monomers. Then a maximum is reached and the chains start to con- 
tract, due to counterion condensation, which decreases the effective charge 
on the chain. When the Bjerrum length is further increased, eventually all 
counterions condense, and then the collapsed chains will start to coagulate 
forming a liquid sol and droplet. This was later also observed in Ref. 31]. 
Because the hydrophobic polyelectrolytes have an additional attractive force 
between the monomers, a reduction of the chain extensions occurs already 
for lower Xb, and this contraction will be stronger than in the hydrophilic 
case. Figure [| reflects exactly this behavior, where we compare the end-to- 
end distances R en d for the two cases e^j = 0.25 and e^j = 0.5 as a function 
of the Manning ratio £m at a density of p = 2 ■ 10 _2 o"~ 3 . The data for the 
hydrophilic polyelectrolytes is taken from Ref. Q. where the model con- 
tained no explicit neutral monomers, so we can not compare the absolute 
values. For our purposes, however, a comparison of the tendencies is suffi- 
cient. In Figure || one recognizes that the hydrophobic system contracts to 
a radius that is a factor of two smaller at £m = 3.94! The relative contrac- 
tion of the hydrophilic chain compared to the reference state is only about 
30% at £m = 7.58. These values clearly show that already the influence 
of weak hydrophobicity changes drastically the conformational properties of 
polyelectrolytes . 

Further evidence for this can be gained from the number of condensed 
counterions as a function of the Bjerrum length Xb, displayed in Figure ^. 
It increases from 12 to 29 around Xb ~ 3 (£m ~ 1) showing exactly the 
dramatic effect as argued above and in reference [fjj]. An analysis of the 
structure factor at the very high Bjerrum length of Xb = 13a, compare 
Figure f|, gives further evidence that the conformation is in a collapsed 
state. 

As seen above, even at a value of eu = 1.0 the chain conformation 
at p = 10 -4 cT -3 was not significantly altered although the neutral chain 
would be strongly in the hydrophobic regime. For that reason, we found 
it necessary to compute an "effective" 9 point for the charged system in 
the following way: The end-to-end distance R en d is calculated as a function 
of The effective 6 point is then defined by the value of e^je ^ at 
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which R en d equals the one for an ideal neutral chain of this length. The 6 
point is called effective, because it is determined only for the system under 
consideration, and hence is not an asymptotic statement. Figure ^ shows the 
time development of R en d for different elj at the fixed density p = 10 _4 cr~ 3 . 

(c) 

From the data the relevant parameter € Ljg ^ is determined to be 

eg,*^ 1.51 ±0.03. (9) 
Polymers with this parameter will be henceforth called "strongly hydropho- 

(c) 

bic". Starting from this value of (- L jq , dominating hydrophobicity at 
higher densities and dominating electrostatics at lower densities will be ob- 
served in the next section. 



4 Variation of the density - single chain properties 



Our intention in this section is to study the crossover from strong screening 
to dominating electrostatic interaction by changing the polymer density p. 
The Bjerrum length is fixed to Xb = 3. Oct, /cgT = 1 and the strength of 
the hydrophobicity is set to e^j = 1.5 in light of the results of the previous 
section. 16 chains of length N = 94 are simulated for densities varying from 
p = 2 ■ KHo-- 3 to p = 2 ■ W- 6 a~ 3 . 

First we want to start by comparing a hydrophilic system to the weakly 
hydrophobic one, e^j = 0.5, and to the strongly hydrophobic system with 
€lj = 1.5. Especially noteworthy is the similarity between the weakly hy- 
drophobic and the hydrophilic system, and the strong contrast of both sys- 
tems to the strongly hydrophobic system. On a very local scale there is no 
large difference in all three systems as a function of density. This is reflected 
in Figure |9| that shows the average bond length. The differences in the total 
values are due to the slight differences in the repulsive part of the LJ po- 
tentials with different clj- As this difference is only small and completely 
irrelevant for the rest of the investigation, we did not correct for it, e.g. by 
an adjustment of the FENE potential. 

Much more interesting is the behavior of the end-to-end distance R en d 

R 2 

and the characteristic ratio r = ^ d . As already presented elsewhere [20] the 
hydrophilic system stretches for decreasing density from a starting point well 
above the self-avoiding-walk value. The hydrophobic chains pass through a 
pronounced minimum at p = 2 • 10~ 2 ct~ 3 . Both observables can be inspected 
as a function of density in Figure 1C. The characteristic ratio r is known to 
equal six for a random walk, 12 for a rod and is supposed to be around 6.3 
for a self avoiding walk. For a homogeneous sphere, with the ends randomly 
distributed over the volume, r = 2, and r = 10/3, if the ends are constraint 
to be randomly distributed on the surface 0, [2(J. For deformed collapsed 
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chains we thus expect a value well below r = 6 but above r = 10/3. For the 
density p = 2 • 10~ 2 o" -3 we find a value of r = 4.68, meaning that the system 
is in a strongly collapsed state. A typical conformation for that density and 



the snapshot of the whole box can be viewed in Figure |1_1. 

At the highest density p = 0.2<t -3 , where the electrostatic interaction is 
even weaker, the chains are not collapsed completely. However, the density 
is so high that strong intermolecular hydrophobic forces cause the system to 
collapse into a gel-like object, as can be inspected in Figure 12, where every 
chain is interpenetrated by several others. The conformations of the single 
chains are dominated by the topological constraints imposed by neighbor- 
ing chains, and we find a broad spectrum of shapes. A simple screening 
argument for the electrostatic and the excluded volume interaction can not 
explain the increase in size, because the density is still too low. However, 
there might also be problems connected with the simulation of a quasi-stable 
(eventually also in experiment "glassy") state. For these higher densities, 
there is almost no energetic advantage to collapse in well separated single 
molecules any more. Figure [13| shows a typical representative of the domi- 
nant conformations. It is clearly less collapsed than the chain at the lower 



density in Figure 11. This can also be observed more clearly in the single 
chain structure factor. A collapsed chain should have a chain length depen- 
dence of R en d that scales with an exponent v = | which leads to a decay of 
the spherically averaged structure factor with a slope of m = —3. This is 
roughly seen in the data set at p = 0.2cr~ 3 . For p = 2 ■ 10 _2 cr _3 , however, 
a value of m = —4 is found that is characteristic for Porod scattering at 
distinct phase boundaries [p^| . This demonstrates how extreme the chain 
collapse is, as can be inspected in Figure [L4|. For lower densities the slope 
converges towards m = —2, which is the random walk value. Nevertheless, 
this value has to be taken with care, because the attractive forces will pre- 
vent one from observing a real random walk behavior. Additionally, the 
linear range in the log-log plot of Figure 14 is rather small. 

We have found that the collapse of the system into a microgel, or iso- 
lated chains, is just a consequence of a strongly hydrophobic attraction in 
polyelectrolytes. Because the simulated densities were rather high, this ef- 
fect should be easily observed in experiments. Especially the isolated chains 
should be a very interesting candidate to study, because they still interact 



strongly via electrostatic interactions. As Figure 11 shows, they should be- 
have like a solution of charge stabilized colloidal particles. Nevertheless, 
the internal degrees of freedom will always play a very important role in 
polyelectrolytes and distinguish them from charged colloids. 

The stability of the colloidal phase can be inferred by the following simple 
estimates p0[|. We start to determine the number of condensed ions in order 
to get the net charge of the collapsed object. The results are plotted in 
Figure 15. We find that 26 ions are condensed at p = 2 ■ 10~ 2 o"~ 3 leading to 
a net charge of 6 per chain. This coincides very nicely with the number which 
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one obtains from the critical charge Q for a charged drop according to the 
Rayleigh criterion, where one balances the Coulomb energy E Cou1 {Rg) = 
X B k B TQ 2 /R G with its surface energy E sur f{R G ) = AnR G k B T/A surf . If 
one uses for the surface area A sur f the monomer area 7r6 2 , one gets for the 
critical charge the formula 

Using the values b = 1.08a and R G = 3.2a, valid for the density p = 2 ■ 
10~ 2 <7~ 3 , one finds Q ~ 6.1. Of course the number of condensed counterions 
strongly increases with density, because the screened electrostatic interaction 
cannot stretch the chains any longer. These more compact states have to 
be stabilized by the oppositely charged counterions. Their entropy loss is 
overcome by the energy gain in the chain. In Ref. [2C] a simple free energy 



argument was presented to show that the colloidal phase (lc) is more stable 
than a hypothetical agglomerate of two chains (2c). The critical charge Qic 
for the agglomerate of two chains for the same density p = 2 • 10 -2 cr~ 3 is 
found with Eq.([To|) by using (R G )2c = Rg^ 1 ^ 3 = 4. Oct to take on the value 
Q2c ~ 8.7. Therefore, in order for two chains to agglomerate an additional 
number of 3.5 counterions have to condense on average onto the pair of 
chains which carry each an average charge of Q c \ = 6.1. To estimate the 
free energy difference of the two globular states, one needs to look only at 
the Coulomb energy of the globules with each other and the entropy of the 
remaining counterions in solution, because the intra Coulomb energy and 
the surface energy are already balanced from the Rayleigh criterion (|Toj) . 
Assuming a uniform distribution of the chains in the simulation box, we 
find the distances d\ c = 18.5a and d.2 C = 23.3cr for the one chain and two 
chain globular state, respectively, to their next neighbors. The repulsive 
Coulomb energy between two neighboring chains is given by 

E Coul = AbA . bT ^_ (n) 

To make the calculation simple we assume the chains to sit on a simple 
cubic lattice. We also neglect higher multi-pole moments and take only 
the next neighbor chains into account. The influence of the counterions 
is also ignored, because on the average only one is available per neighbor. 
With these simplifications the average Coulomb interaction energy per single 
chain is found to be E*(° ul = 18.1k B T '/chain and E%° ul = 14.6fc B T/chain, 
respectively. The corresponding value for an agglomerate of two chains is 
therefore about 3.5k B T lower, preferring the agglomerate. The entropy of 
the free counterions can be approximated by an ideal gas Ansatz. The free 
counterion densities are given by p\ c = 6.1/(94 + 32)p and p2 C = 4.35/ (94 + 
32) p, respectively. Because the entropy per free ion is given by the logarithm 
of the available volume fraction the entropic contribution of the free energy 
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difference per single chain is given by 6.11npi c — 4.351np2c ~ — 10.7/sbT. 
This gives a total preference of 7. 2k bT /chain for the single chain globule |2C]. 

The spherically averaged single chain structure factor shows another very 
striking feature at low q, this time for the lowest densities. At p = 2-10~ 5 cr~ 3 
and p = 2 • 10~ 6 <t~ 3 a second well defined length scale shows up. The log- 
arithmic slope measured is exactly m = — 1, as can be seen in Figure 
which corresponds to an exponent u = 1. This is an extraordinary result, 
because no other polymeric systems show such a behavior. For example, hy- 
drophilic polyelectrolytes always have an exponent v < 1, because entropic 
fluctuations lead to some local roughness 0. Even the blob pole conforma- 
tions observed in simulations of single Debye-Hiickel chains Q that show 
a linear dependency of the end-to-end distance on the contour length after 
introduction of the concept of electrostatic blobs, show an effective expo- 
nent v smaller than one. Figure 17 and 18 show typical pearl-necklace like 
conformations. Small globules are connected by thin bridges. The stronger 
the electrostatic interactions are, the larger is the number of globules. The 
globules are in addition stabilized by a few condensed counterions. So we 
can qualitatively reproduce the picture of Ref. [JDJ for our strongly charged 
systems. Of course, the globules are much smaller than in their weakly 
charged case, but clearly detectable. The chains are strongly stretched due 
to the electrostatic interaction which orients the different globules to an 
almost rodlike object by maximizing their mutual distance. Additionally, 
the chains are, at least within the globules, locally not very flexible, be- 
cause the monomers dislike to break the energetically favorable hydrophobic 
bonds (around E = 1.5/jbT per contact). To complete the set of pictures 
of conformations, we display in Figures [l^ and |2(] typical conformations at 
intermediate densities. A continuous deformation towards a dumbbell like 
conformation is observed, but no cylindrical instability is found. It is inter- 
esting to note, that the larger globules are found at the chain ends, and the 
smaller ones reside in the chain middle. 

At the end of this section we want to stress again, that no micro-phase 
separation is observed, although this is theoretically predicted (JOJ. Only 
for the highest density we can not rule out that this will eventually happen. 
The main argument for a micro-phase separation was that the entropy of 
the condensed counterions in the polyelectrolyte-rich phase is much larger, 
because the condensed counterions can explore the whole volume of the 
separated phase. We have been showing in our previous communication 
[20 1, that the system prefers single chain globules, which repel each other 
because they still have a sufficiently high net charge. Additional counterions 
would need to condense to stabilize an agglomerate of two or more chains. 
However, their loss of translational entropy (and the increase in Coulomb 
energy) is much larger than the gain in surface energy, and therefore single, 
isolated globules are clearly favored in certain density regimes. 



11 



5 Comparison with experiments 



It is important to try to compare these results to experimental data. Un- 
fortunately, most experiments do not examine single chain properties with 
adequate precision. Especially for low densities, where the mutual influence 
of the chains can be neglected, scattering experiments are extremely diffi- 
cult. A comparably new experimental study [16| shows good agreement with 
our data at relatively high densities. Unfortunately, the density was only 
varied once by a factor of two, so that the effects we report in this paper, 
could not be seen. 

Most experimental results come from detailed measurements of the struc- 
ture factor of the whole sample or thermodynamic properties like the osmotic 
pressure II are measured. To compare our results of section^, however, data 
of single chain properties are needed, which unfortunately are rare. The 
structure factor of the whole simulation box that reflects the arrangement 
of the chains in space, is difficult to determine from the simulation data 
for the present poor solvent situation. The reasons are of technical nature. 
The structure factor itself can easily be calculated, however it is strongly 
influenced by finite size effects, as the number of chains in the simulation 
box is relatively small (16 or 32). The collapsed chains can see each other 
only weakly, which results in broad maxima in the pair-correlation function 
and in S(q), which was not the case for the hydrophilic systems Nev- 
ertheless, the simulations are able to show some important aspects, which 
can help to interpret the experimental results, as will be shown below. 

Experimentally, the maximum of the structure factor, which is the Fourier 
transform of the monomer pair correlation function of all monomers in the 
sample, shows a concentration dependent position q max ~ Pc f with v s f = ^ 
at high densities and v s f = | at low densities, where p c denotes the charge 
concentration. Because the statistical fluctuations are too high, a reason- 
able distinction between the relevant exponents (0.5 to 0.3) can not be given 
in our case. At low densities the scattering should reflect nothing but the 
distances between the uniformly distributed chains. These distances vary 
with p^ 1 / 3 . At higher densities the scattering results from the correlation 
hole that is supposed to be proportional to the screening length which then 
is proportional to the square root of the charge concentration J?], [l5|, 16 1. 

For hydrophobic polyelectrolytes deviations of this behavior are found in 
small angle neutron and x-ray scattering experiments fi~6| , 33|. An effective 
exponent of u s f = 0.4 is observed, while the hydrophilic chains show the 
usual exponent v s f = 0.5. The investigated hydrophilic polymer is strongly 
charged and seems to follows the Manning criterion over the examined small 
density range. So the charge is "fixed" , meaning the number of free charges 
in the system remains constant. In the following we argue that the num- 
ber of free charges in the system is the relevant parameter to explain the 
experimental results. In Figure 15 the number of condensed counterions 
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was plotted, which showed a strong density dependence. Because every 
condensed counterion neutralizes effectively one charged monomer, this will 
lead naturally to a reduced charged density at higher mass densities. The 
calculation of the screening length based on this "effective" charge density 
reproduces the smaller exponents quite well. Only at the highest densities, 
where the number of condensed counterions is almost constant, the exponent 
is given by the volume change and equals u s f = 0.5. In Table |l| and [2] we 
estimate the effective hard sphere diameter of the chains by calculating the 
DH radius. We assume that screening is only facilitated by the mobile, non- 
condensed counterions with a density pt c . This defines a screening length 
= (AirXsPfc) 1 ^ 2 ■ If w e assume that the solution is in a liquid of hard 
spheres of radius k™ 1 and if we further assume that the effective particles 
are certainly big enough to interact, than we can argue that the effective 
exponent y, obtained via oc p y as given in Tables || and || should reflect 
the experimental data for the scattering. There is however a word of care 
needed. The value of v s f = 0.5 as given from scaling assumes a blob picture 
with strongly interacting chains. Thus this cannot expected to hold so easily 
here. On the other hand for hydrophilic, shorter chains the crossover from 
v s f = 1/3 to u s f = 1/2 was observed [f7|]. 

Unfortunately quite often one considers in the experimental literature 
only the density of the chemical groups that can dissociate a counterion. 
We term this quantity the "chemical" charge density in this section, because 
it is the maximal charge density possible for chemical reasons. The actual 
charge density will almost always be smaller due to Manning condensation 
and additionally condensation due to collapse. 

In Ref. | |15| |, for example, a charge density of / = 0.04 was found for a 
NaPSS molecule by osmotic pressure measurements, with a chemical charge 
density given by f c hem = 0.38. Although all analyzed systems are well above 
the critical Manning value of one, and should therefore show a fixed charge 
density, a strong variation of the measured charge densities was observed. 
This is a clear hint that strong condensation effects due to more compact 
conformations play an important role for hydrophobic polyelectrolytes, in 
accordance to the simulation results of section ||[ This has also consequences 
for the osmotic pressure [34]. 

The determination of the fraction of bound (immobile) counterions is 
very difficult to accomplish experimentally. We imagine that particularly 
with x-ray and synchrotron scattering for some well prepared systems one 
could achieve a sufficient contrast between the counterions and the chains. 
As the majority of the chain atoms are carbon, hydrogen and oxygen the use 
of heavy counterions like cesium should give a sufficient contrast in scatter- 
ing intensities. Therefore the structure factor of the counterion distribution 
should be measurable. This quantity is much easier calculated in the simu- 
lation, because the counterions are more mobile and outnumber the number 
of polymers by far. So their distribution in the system can be dissolved 
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in much more detail. The simulations show very interesting results, as for 
example in Figure 21 where the density dependency of the small q part of 
the counterion structure factor is shown. For high and medium densities 
the distribution is coupled to that of the polymers and shows a pronounced 
signal of the underlying structure. For the highest density, the structure 
factor shows no clear peak. This is the reflection of the relatively uniform 
distribution of the counterions in the microgel. In the medium density range 
a maximum shows up which increases and moves to smaller g-values with 
decreasing density. In the case of p = 2 ■ 10~ 5 o"~ 3 the intensity decreases, 
because more and more counterions leave the strongly stretched chains so 
that a significant background contributes to the scattering. For the lowest 
density only this background of uniformly distributed counterions remains, 
leading to a completely flat and unstructured scattering function. For higher 
q a clear minimum is detected, which is getting deeper with decreasing den- 
sity, as can be seen in Figure Of course, the lowest density still gives the 
flat signal from the uniform background. This minimum can be interpreted 
as a correlation hole between the counterion clouds. At higher densities, 
many of the counterions are condensed and therefore near the chains. The 
number of free counterions is relatively small compared to the free volume 
between the chains which spans almost the whole volume of the solution, 
as can be seen in Figure 11. This concentration difference is reflected by 
the structure factor. At even higher q- values a second maximum can be 
detected for the higher densities in the simulation data. It is located around 
the g-value related to the radius of gyration of the polyelectrolytes reflect- 
ing the correlation of the chains within the counterion cloud around one 
monomer. As the amplitude is very small this effect is difficult to analyze 
quantitatively, however we consider it not to be relevant for experimental 
studies at present, and do not analyze it in further detail here. 



6 Conclusions 

We present a set of MD simulations studying the influence of a poor sol- 
vent on the structural properties of strongly charged, flexible polyelectrolyte 
chains in salt free solution. 

We showed that the underlying neutral system is characterized by a 6 
point which we determined to be ehjfi = 0.34 ± 0.02. In the beginning 
three different values of the strength of the hydrophobicity (eu = 0.25, 
0.5, and 1.0) were examined at a polymer density of p = 10 -4 <7~ 3 . The 
value of €lj = 0.25 lead to hydrophilic behavior as expected, because this 
value is even below the neutral 9 point. However, the other two values 
of €lj showed as well only a weak influence of the attractive interaction 
and the differences to the purely hydrophilic chains was rather small. By 
varying the Bjerrum length Xb however, we could demonstrate that under 
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the right conditions weak hydrophobicities can nevertheless lead to strong 
effects. For e^j = 0.5 the chains showed a stretching followed by a strong 
contraction with increasing Bjerrum length. This effect is well known from 
simulations of hydrophilic polyelectrolytes, but in our hydrophobic case it 
is more pronounced and shifted to smaller A#- 

An "effective" 9 point was calculated for the charged system for a density 
of p = 10 -4 ct -3 . The condition that its end-to-end distance equals the one 

(c) 

for a ideal, neutral chain of the same length resulted in € Lje ^ = 1.51±0.03. 
This value was taken as a starting point for detailed study of the chain con- 
formations under a variation of the polymer density. At high densities we 
found a high screening of the electrostatic interactions, favoring the short 
range hydrophobic attraction, whereas at low densities the Coulomb forces 
clearly dominated. A multitude of chain conformations ranging from micro- 
gels via isolated, strongly collapsed chains to extremely stretched necklace 
chains was found. The collapsed objects showed locally Porod scattering 
(S(q) oc q~ 4 ), while the necklace like objects showed a logarithmic slope of 
their scaling function of m = — 1 {y = 1), so that both extremes are realized 
in one system that differs only by density. A micro-phase separation was not 
observed. We estimated from our simulation that an agglomeration of two 
chains is unfavorable for the density p = 2 • 10~ 2 cr~ 3 , whereas more detailed 
estimates are given in |2C]. To sum up, the polymer density seems to be a 
very useful parameter to vary to study polyelectrolytes under poor solvent 
conditions. The preparation of such systems should be extremely fruitful 
also for experiments, because the density p can be very easily controlled. 
The comparison of the MD data of the structure factor of the monomer pair 
correlation with experimental data indicates that the knowledge of the real 
number of free charges in the system is of central importance in understand- 
ing the behavior of polyelectrolyte solutions, especially for the poor solvent 
case. Measurements of the counterion structure factor would be very helpful 
in comparison to the simulation results presented here. 

We could relate the structural properties of polyelectrolytes under good 
and poor solvent conditions within one consistent approach. This gave a 
clear and conclusive picture of the influence of solvent quality on the struc- 
ture of flexible polyelectrolyte chains. 
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List of Abbrevations and Symbols 

MD Molecular Dynamics 

LJ Lennard-Jones 

FENE finite extension nonlinear elastic 

PME Particle-Mesh-Ewald 

PSS sulfonated polystyrene 

U L j,hydrophil purely repulsive hydrophilic potential 

ULJ,hydrophob hydrophobic potential showing a short range attraction 

Ufene Finite Extension Nonlinear Elastic bond potential between chain 
monomers 

R en d end-to-end distance 

Rg radius of gyration 

R 2 

r characteristic ratio r = -g^ 

N Numbers of monomers of a chain 

L linear length of the simulation cell 

b average bond length 

Xb Bjerrum length 

k Debye-Hiickel parameter 

e unit charge 

q valence of a unit charge 

/ fraction of charged monomers 

fchem chemical charge density 

Q critical charge determined by the Rayleigh criterion 

E Coui Coulomb energy 

d distance between globules 

ks Boltzmann constant 

T Temperature 

a Lennard-Jones length unit (monomer radius) 

r Lennard-Jones time unit 

eo, es permittivity of the vacuum and the solvent, respectively 

r friction coefficient of the Langevin dynamics 

At time step of the velocity Verlet MD algorithm 

£m Manning Parameter 
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£^ critical Manning parameter 

p density: number of monomers per simulation volume 
9 - temperature of the neutral system 

(c) 

Tq ' 9 - temperature of the charged system 

(n) 

€ LJ6 energy parameter giving the 9 point of the neutral system 

£ lj e f f energy parameter giving the effective 9 - point of the charged 
system 

T reduced temperature to the critical 9 - point 

S sp sphericall averaged structure factor 

S\\ component of the spherically averaged structure factor along 
the first principal axis of inertia 

q magnitude of the Fourier space vector 

Qmax position of the maximum of the structure factor 

II osmotic pressure 

v scaling exponent of the chain extension 

v s f scaling exponent of the maximum of the structure factor 

p c charge concentration 

pf c density of mobile, noncondensed counterions 
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Variation of the effective screening for hydrophobic polyelectrolytes 
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Table 1: Variation of the effective screening with density for hydrophobic 
polyelectrolytes. Nc denotes the number of condensed counterions, N c h is 
the number of free charges, and p c h is the density of free charges. The 
exponent y is obtained via oc p y . One recognizes a deviation of y = 1/2 
due to counterion condensation 



Variation of the effective screening for hydrophilic polyelectrolytes 
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Table 2: Variation of the effective screening with density for hydrophilic 
polyelectrolytes. The notation is the same as in Table ||. The variations are 
clearly much weaker than for the hydrophobic case. 
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Figure 1: Scaling analysis for determining the point of the neutral system. 
Plotted is e^j vs. N 05 according to Eq. [/], where T — e L j in our units. 

One obtains straight lines N~2 — > which converge up to corrections of 
higher order to the point. 6 e ff describes the point for the finite chain, 
where e^j was chosen such that its R en d had the same value as an equally 
long random walk. These values converge as well to the 9 point for N — > oo. 
One obtains t£j e = 2 96 xg 6 = ^ 0-02. 
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Figure 2: Comparison of the spherically averaged structure factor, S sp (q) 
(top curve) and its component along the first principal axis of inertia, S\\(q) 
(bottom curve), versus q for various hydrophobic interaction strengths e^j = 
0.25, 0.5, and 1.0 at a density p = 10-V" 3 . 
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Figure 3: Comparison of the radius of gyration Rq versus density p of 
hydrophilic (phil) and weakly hydrophobic (elj = 0.5) polyelectrolytes. 
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Figure 4: The spherically averaged structure factor S sp (q) vs. q for hy- 
drophilic (phil) and weak hydrophobic (e^j = 0.5) chains for the density 
p = 10 _1 (T -3 . The dashed straight lines indicate the slope of the linear part 
of the hydrophobic data, the solid lines that one of the hydrophilic data. 
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Figure 5: End-to-end distance R en d as a function of the Manning parameter 
£m for hydrophilic chains (eu = 0.25, □, from Ref. Q), and hydrophobic 
chains (e^j = 0.5, •). The straight line indicates comparable end-to-end 
distances. 
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Figure 6: Number of condensed counter ions Nc as a function of the Bjerrum 
length \b for eij = 0.5 at the density p = 2 • lCT 3 cr~ 3 . The number of 
counterions per chain is 32. 
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Figure 7: Spherically averaged structure factor S sp (q) vs. q for various 
Bjerrum lengths Xb for hydrophobic chains with elj = 0.5. 
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Figure 8: Time development in MD steps of the end-to-end distance R en d 
of charged chains for various hydrophobicities elj- The straight horizontal 
line indicates the end-to-end distance of the equally long uncharged random 
walk 
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Figure 9: Bond length b versus density p for hydrophilic (phil), weakly hy- 
drophobic (eu = 0.5), and hydrophobic (eij = 1.5) polyelectrolyte chains. 
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Figure 10: Density dependence of the end-to-end distance R en d (top) and 
the statistical ratio r (bottom) vs. density p for hydrophilic (phil), weakly 
hydrophobic (e^j = 0.5), and strongly hydrophobic (clj = 1.5) systems. 
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Figure 11: Conformations at the density p = 2 • 10~ 2 cr~ 3 . Top: A typi- 
cal poor solvent polyelectrolyte conformation. The beads are the charged 
monomers, the counterions are indicated as small black spheres, and the 
neutral monomers are dark grey spheres. Only counterions within a dis- 
tance of 3<7 to the chain are displayed. Bottom: A snapshot of the whole 
simulation box, showing all 16 chains together with their counterions. The 
picture shows that the chains collapse into single globules which are well 
separated, and a small fraction of the counterions is still in solution. 
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Figure 12: A snapshot of the whole simulation box, showing 16 chains to- 
gether with their counterions for the density p = 2 • 10 _1 o" -3 . The beads are 
the charged monomers, the counterions are indicated as small black spheres, 
and the neutral monomers are dark grey spheres. 
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Figure 14: Details of the spherically averaged structure factor S sp (q) for large 
q for the indicated densities p. The straight lines have slope m = —2, —3, 
and —4, respectively. 
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Figure 15: The number N c of condensed counter ions per chain as a function 
of density p for 32 charged monomers per chain for hydrophilic chains (phil) 
and hydrophobic chains with eij = 0.5 and 1.5. 
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Figure 16: Details of the spherically averaged structure factor S sp (q) vs. q 
for small q, obtained for the indicated densities p. The straight lines have 
all slope m = — 1, and the vertical lines denote the value of q = 2-k/Rg- 



36 




Figure 17: A typical poor solvent polyelectrolyte conformation for the den- 
sity p = 2 ■ lCP 5 cr~ 3 , which shows nicely a pearl-necklace structure. The 
coloring is as in Fig. 11 
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Figure 19: A typical poor solvent polyelectrolyte conformation for the den- 
sity p = 2 • 10~ 3 a~ 3 . One sees the onset of the Rayleigh instability, the 



globule wants to split up into two parts. The coloring is as in Fig. 11 
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Figure 20: A typical poor solvent polyelectrolyte conformation for the den- 
sity p = 2 ■ 10 -4 cr -3 , where now two globules with a bridge connecting them, 
has formed . The coloring is as in Fig. 11 
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Figure 21: The structure factor Si on (q) vs. q of the counter ion distribution 
for the indicated densities p. For lower densities the structure factor becomes 
similar to the one of a homogeneous background. 
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Figure 22: Blow up of Fig. 21 for large q. One can detect a minimum of the 
structure factor for q ~ 1.0, which is reminiscent of a correlation hole of the 
counterion cloud. 
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